import numpy as np




def source_vector(vspace,PDE):
     mesh = vspace.mesh
     pp = mesh.entity_barycenter('cell')
     
     alpha = PDE.get_alpha(pp)
     idx, = np.nonzero((alpha < 2*np.pi/9)&(alpha>0.0))
     

     cellmeasure = vspace.cellmeasure
     bcs, ws = vspace.integrator.get_quadrature_points_and_weights()
     rho = mesh.bc_to_point(bcs)[:,idx,0] #(NQ,NC)

     gdof = vspace.number_of_global_dofs()
     shape = (gdof,2)
     b = np.zeros(shape,dtype=np.float64) #(gdof,2)
     phi = vspace.basis(bcs) #(NQ,NC,ldof)



     bb = np.einsum('m,mi,mik,i->ik',ws,rho,phi,cellmeasure[idx])*PDE.F #(NC0,ldof)
     cell2dof = vspace.cell_to_dof()[idx] #(NC0,ldof)


     np.add.at(b[:,1], cell2dof, bb)
    

     return b
